library(sf)
library(ggplot2)
library(leaflet)
library(scales)
library(ggmap)
library(dplyr)
ak_regions <- read_sf("../shapefiles/ak_regions_simp.shp")
plot(ak_regions)
class(ak_regions)
## [1] "sf" "tbl_df" "tbl" "data.frame"
head(ak_regions)
## Simple feature collection with 6 features and 3 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: -179.2296 ymin: 51.15702 xmax: 179.8567 ymax: 71.43957
## Geodetic CRS: WGS 84
## # A tibble: 6 x 4
## region_id region mgmt_area geometry
## <int> <chr> <dbl> <MULTIPOLYGON [°]>
## 1 1 Aleutian Islands 3 (((-171.1345 52.44974, -171.1686 52.4174~
## 2 2 Arctic 4 (((-139.9552 68.70597, -139.9893 68.7051~
## 3 3 Bristol Bay 3 (((-159.8745 58.62778, -159.8654 58.6137~
## 4 4 Chignik 3 (((-155.8282 55.84638, -155.8049 55.8655~
## 5 5 Copper River 2 (((-143.8874 59.93931, -143.9165 59.9403~
## 6 6 Kodiak 3 (((-151.9997 58.83077, -152.0358 58.8271~
ak_regions_3338 <- ak_regions %>% #use pipe because SF is part of tidyverse/dplyr
st_transform(crs = 3338)
st_crs(ak_regions_3338)#this is a good CRS projection for AK
## Coordinate Reference System:
## User input: EPSG:3338
## wkt:
## PROJCRS["NAD83 / Alaska Albers",
## BASEGEOGCRS["NAD83",
## DATUM["North American Datum 1983",
## ELLIPSOID["GRS 1980",6378137,298.257222101,
## LENGTHUNIT["metre",1]]],
## PRIMEM["Greenwich",0,
## ANGLEUNIT["degree",0.0174532925199433]],
## ID["EPSG",4269]],
## CONVERSION["Alaska Albers (meters)",
## METHOD["Albers Equal Area",
## ID["EPSG",9822]],
## PARAMETER["Latitude of false origin",50,
## ANGLEUNIT["degree",0.0174532925199433],
## ID["EPSG",8821]],
## PARAMETER["Longitude of false origin",-154,
## ANGLEUNIT["degree",0.0174532925199433],
## ID["EPSG",8822]],
## PARAMETER["Latitude of 1st standard parallel",55,
## ANGLEUNIT["degree",0.0174532925199433],
## ID["EPSG",8823]],
## PARAMETER["Latitude of 2nd standard parallel",65,
## ANGLEUNIT["degree",0.0174532925199433],
## ID["EPSG",8824]],
## PARAMETER["Easting at false origin",0,
## LENGTHUNIT["metre",1],
## ID["EPSG",8826]],
## PARAMETER["Northing at false origin",0,
## LENGTHUNIT["metre",1],
## ID["EPSG",8827]]],
## CS[Cartesian,2],
## AXIS["easting (X)",east,
## ORDER[1],
## LENGTHUNIT["metre",1]],
## AXIS["northing (Y)",north,
## ORDER[2],
## LENGTHUNIT["metre",1]],
## USAGE[
## SCOPE["Topographic mapping (small scale)."],
## AREA["United States (USA) - Alaska."],
## BBOX[51.3,172.42,71.4,-129.99]],
## ID["EPSG",3338]]
plot(ak_regions_3338)
ak_regions_3338 %>%
select(region)
## Simple feature collection with 13 features and 1 field
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: -2175328 ymin: 405653 xmax: 1579226 ymax: 2383770
## Projected CRS: NAD83 / Alaska Albers
## # A tibble: 13 x 2
## region geometry
## <chr> <MULTIPOLYGON [m]>
## 1 Aleutian Islands (((-1156666 420855.1, -1159837 417990.3, -1161898 41694~
## 2 Arctic (((571289.9 2143072, 569941.5 2142691, 569158.2 2142146~
## 3 Bristol Bay (((-339688.6 973904.9, -339302 972297.3, -339229.2 9710~
## 4 Chignik (((-114381.9 649966.8, -112866.8 652065.8, -108836.8 65~
## 5 Copper River (((561012 1148301, 559393.7 1148169, 557797.7 1148492, ~
## 6 Kodiak (((115112.5 983293, 113051.3 982825.9, 110801.3 983211.~
## 7 Kotzebue (((-678815.3 1819519, -677555.2 1820698, -675557.8 1821~
## 8 Kuskokwim (((-1030125 1281198, -1029858 1282333, -1028980 1284032~
## 9 Cook Inlet (((35214.98 1002457, 36660.3 1002038, 36953.11 1001186,~
## 10 Norton Sound (((-848357 1636692, -846510 1635203, -840513.7 1632225,~
## 11 Prince William Sound (((426007.1 1087250, 426562.5 1088591, 427711.6 1089991~
## 12 Southeast (((1287777 744574.1, 1290183 745970.8, 1292940 746262.7~
## 13 Yukon (((-375318 1473998, -373723.9 1473487, -373064.8 147393~
pop <- read.csv("../shapefiles/alaska_population.csv")
head(pop)
## year city lat lng population
## 1 2015 Adak 51.88000 -176.6581 122
## 2 2015 Akhiok 56.94556 -154.1703 84
## 3 2015 Akiachak 60.90944 -161.4314 562
## 4 2015 Akiak 60.91222 -161.2139 399
## 5 2015 Akutan 54.13556 -165.7731 899
## 6 2015 Alakanuk 62.68889 -164.6153 777
class(pop)
## [1] "data.frame"
- Do a spatial join using st_join() - look at documentation for st_join for other joins
pop_4326 <- st_as_sf(pop,
coords = c('lng', 'lat'),
crs = 4326, #4236 is WGS84. You have to read in the data in the format it's in and then transform it. You cannot just assign it a diff CRS (such as 3338)
remove = F)
head(pop_4326)
## Simple feature collection with 6 features and 5 fields
## Geometry type: POINT
## Dimension: XY
## Bounding box: xmin: -176.6581 ymin: 51.88 xmax: -154.1703 ymax: 62.68889
## Geodetic CRS: WGS 84
## year city lat lng population geometry
## 1 2015 Adak 51.88000 -176.6581 122 POINT (-176.6581 51.88)
## 2 2015 Akhiok 56.94556 -154.1703 84 POINT (-154.1703 56.94556)
## 3 2015 Akiachak 60.90944 -161.4314 562 POINT (-161.4314 60.90944)
## 4 2015 Akiak 60.91222 -161.2139 399 POINT (-161.2139 60.91222)
## 5 2015 Akutan 54.13556 -165.7731 899 POINT (-165.7731 54.13556)
## 6 2015 Alakanuk 62.68889 -164.6153 777 POINT (-164.6153 62.68889)
class(pop_4326)
## [1] "sf" "data.frame"
- Oops an error! Because they use different projections
#pop_joined <- st_join(pop_4326, ak_regions_3338, join = st_within) comment out bc not working, don't need?
pop_3338 <- st_transform(pop_4326, crs = 3338)
head(pop_3338)
## Simple feature collection with 6 features and 5 fields
## Geometry type: POINT
## Dimension: XY
## Bounding box: xmin: -1537925 ymin: 472626.9 xmax: -10340.71 ymax: 1456223
## Projected CRS: NAD83 / Alaska Albers
## year city lat lng population geometry
## 1 2015 Adak 51.88000 -176.6581 122 POINT (-1537925 472626.9)
## 2 2015 Akhiok 56.94556 -154.1703 84 POINT (-10340.71 770998.4)
## 3 2015 Akiachak 60.90944 -161.4314 562 POINT (-400885.5 1236460)
## 4 2015 Akiak 60.91222 -161.2139 399 POINT (-389165.7 1235475)
## 5 2015 Akutan 54.13556 -165.7731 899 POINT (-766425.7 526057.8)
## 6 2015 Alakanuk 62.68889 -164.6153 777 POINT (-539724.9 1456223)
pop_joined <- st_join(pop_3338, ak_regions_3338, join = st_within)
head(pop_joined)
## Simple feature collection with 6 features and 8 fields
## Geometry type: POINT
## Dimension: XY
## Bounding box: xmin: -1537925 ymin: 472626.9 xmax: -10340.71 ymax: 1456223
## Projected CRS: NAD83 / Alaska Albers
## year city lat lng population region_id region
## 1 2015 Adak 51.88000 -176.6581 122 1 Aleutian Islands
## 2 2015 Akhiok 56.94556 -154.1703 84 6 Kodiak
## 3 2015 Akiachak 60.90944 -161.4314 562 8 Kuskokwim
## 4 2015 Akiak 60.91222 -161.2139 399 8 Kuskokwim
## 5 2015 Akutan 54.13556 -165.7731 899 1 Aleutian Islands
## 6 2015 Alakanuk 62.68889 -164.6153 777 13 Yukon
## mgmt_area geometry
## 1 3 POINT (-1537925 472626.9)
## 2 3 POINT (-10340.71 770998.4)
## 3 4 POINT (-400885.5 1236460)
## 4 4 POINT (-389165.7 1235475)
## 5 3 POINT (-766425.7 526057.8)
## 6 4 POINT (-539724.9 1456223)
- we can drop geometries now, since we no longer need it. But since GEOMETRIES are sticky you need to unstick it using “as.data.frame” to convert if from SF table to normal dataframe
pop_region <- pop_joined %>%
as.data.frame() %>%
group_by(region) %>%
summarise(total_pop = sum(population))
head(pop_region)
## # A tibble: 6 x 2
## region total_pop
## <chr> <int>
## 1 Aleutian Islands 8840
## 2 Arctic 8419
## 3 Bristol Bay 6947
## 4 Chignik 311
## 5 Cook Inlet 408254
## 6 Copper River 2294
- Do a regular left_join to add the total_pop to the original map talbe
pop_region_3338 <- left_join(ak_regions_3338, pop_region)
## Joining, by = "region"
#joining by region
plot(pop_region_3338)
pop_mgmt_3338 <- pop_region_3338 %>%
group_by(mgmt_area) %>%
summarize(total_pop = sum(total_pop), do_union = F) #do union to get rid of holes???
head(pop_mgmt_3338)
## Simple feature collection with 4 features and 2 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: -2175328 ymin: 405653 xmax: 1579226 ymax: 2383770
## Projected CRS: NAD83 / Alaska Albers
## # A tibble: 4 x 3
## mgmt_area total_pop geometry
## <dbl> <int> <MULTIPOLYGON [m]>
## 1 1 67252 (((1287777 744574.1, 1290183 745970.8, 1292940 746262.7, ~
## 2 2 417595 (((561012 1148301, 559393.7 1148169, 557797.7 1148492, 55~
## 3 3 24224 (((-1156666 420855.1, -1159837 417990.3, -1161898 416944.~
## 4 4 137476 (((571289.9 2143072, 569941.5 2142691, 569158.2 2142146, ~
plot(pop_mgmt_3338["mgmt_area"])
- map with ggplot. include args explicitly instead of shorthand, bc geom_sf puts data and mapping args in reverse order!
ggplot(pop_region_3338) +
geom_sf(aes(fill = total_pop))
ggplot(pop_region_3338) +
geom_sf(aes(fill = total_pop)) +
theme_bw() +
labs(fill = "Total Population") +
scale_fill_continuous(low = "khaki", high = "firebrick", labels = comma)
- Need to add point and river data
rivers_3338 <- read_sf("../shapefiles/ak_rivers_simp.shp")
st_crs(rivers_3338)
## Coordinate Reference System:
## User input: Albers
## wkt:
## PROJCRS["Albers",
## BASEGEOGCRS["GCS_GRS 1980(IUGG, 1980)",
## DATUM["D_unknown",
## ELLIPSOID["GRS80",6378137,298.257222101,
## LENGTHUNIT["metre",1,
## ID["EPSG",9001]]]],
## PRIMEM["Greenwich",0,
## ANGLEUNIT["Degree",0.0174532925199433]]],
## CONVERSION["unnamed",
## METHOD["Albers Equal Area",
## ID["EPSG",9822]],
## PARAMETER["Latitude of false origin",50,
## ANGLEUNIT["Degree",0.0174532925199433],
## ID["EPSG",8821]],
## PARAMETER["Longitude of false origin",-154,
## ANGLEUNIT["Degree",0.0174532925199433],
## ID["EPSG",8822]],
## PARAMETER["Latitude of 1st standard parallel",55,
## ANGLEUNIT["Degree",0.0174532925199433],
## ID["EPSG",8823]],
## PARAMETER["Latitude of 2nd standard parallel",65,
## ANGLEUNIT["Degree",0.0174532925199433],
## ID["EPSG",8824]],
## PARAMETER["Easting at false origin",0,
## LENGTHUNIT["metre",1],
## ID["EPSG",8826]],
## PARAMETER["Northing at false origin",0,
## LENGTHUNIT["metre",1],
## ID["EPSG",8827]]],
## CS[Cartesian,2],
## AXIS["(E)",east,
## ORDER[1],
## LENGTHUNIT["metre",1,
## ID["EPSG",9001]]],
## AXIS["(N)",north,
## ORDER[2],
## LENGTHUNIT["metre",1,
## ID["EPSG",9001]]]]
head(rivers_3338)
## Simple feature collection with 6 features and 3 fields
## Geometry type: MULTILINESTRING
## Dimension: XY
## Bounding box: xmin: -622169.8 ymin: 557375.7 xmax: 738652.5 ymax: 2316718
## Projected CRS: Albers
## # A tibble: 6 x 4
## StrOrder region n geometry
## <int> <chr> <int> <MULTILINESTRING [m]>
## 1 3 Aleutian Islands 316 ((-616509.9 557976.9, -616240.8 557375.7), (-~
## 2 3 Arctic 10869 ((-20959.19 2009003, -20804.35 2009054), (-20~
## 3 3 Bristol Bay 4884 ((-231927.5 794449.9, -232004.6 794730.2), (-~
## 4 3 Chignik 141 ((-302874.8 715936.3, -302829.4 715920.8), (-~
## 5 3 Cook Inlet 3447 ((-11441.06 987666.2, -11772.61 988142.6), (-~
## 6 3 Copper River 1784 ((540428.6 1180052, 539410.1 1176802), (54043~
ggplot(pop_region_3338) +
geom_sf(aes(fill = total_pop)) +
geom_sf(data = rivers_3338, mapping = aes(size = StrOrder)) +
theme_bw() +
labs(fill = "Total Population") +
scale_fill_continuous(low = "khaki", high = "firebrick", labels = comma)
- need to scale size of rivers to something skinner - guess and check!
ggplot(pop_region_3338) +
geom_sf(aes(fill = total_pop)) +
geom_sf(data = rivers_3338, mapping = aes(size = StrOrder)) +
geom_sf(data = pop_3338, mapping = aes(), size = 0.5) +
scale_size(range = c(0.01, 0.2), guide = "none") +
theme_bw() +
labs(fill = "Total Population") +
scale_fill_continuous(low = "khaki", high = "firebrick", labels = comma)
- package ggmap let’s you plot your own points on a basemap (base tiles)
- unlike leaflet you can’t zoom around on it.
- basemap is in CRS 3857, so we need to transform our data to match
pop_3857 <- pop_3338 %>%
st_transform(crs = 3857)
- get the basemap, but first fix an issue with ggmap with this function
# Define a function to fix the bbox to be in EPSG:3857
# See https://github.com/dkahle/ggmap/issues/160#issuecomment-397055208
ggmap_bbox_to_3857 <- function(map) {
if (!inherits(map, "ggmap")) stop("map must be a ggmap object")
# Extract the bounding box (in lat/lon) from the ggmap to a numeric vector,
# and set the names to what sf::st_bbox expects:
map_bbox <- setNames(unlist(attr(map, "bb")),
c("ymin", "xmin", "ymax", "xmax"))
# Coonvert the bbox to an sf polygon, transform it to 3857,
# and convert back to a bbox (convoluted, but it works)
bbox_3857 <- st_bbox(st_transform(st_as_sfc(st_bbox(map_bbox, crs = 4326)), 3857))
# Overwrite the bbox of the ggmap object with the transformed coordinates
attr(map, "bb")$ll.lat <- bbox_3857["ymin"]
attr(map, "bb")$ll.lon <- bbox_3857["xmin"]
attr(map, "bb")$ur.lat <- bbox_3857["ymax"]
attr(map, "bb")$ur.lon <- bbox_3857["xmax"]
map
}
- basemap we are using is stamenmap
bbox <- c(-170, 52, -130, 68) # This is roughly southern Alaska
ak_map <- get_stamenmap(bbox, zoom = 4)
## Source : http://tile.stamen.com/terrain/4/0/3.png
## Source : http://tile.stamen.com/terrain/4/1/3.png
## Source : http://tile.stamen.com/terrain/4/2/3.png
## Source : http://tile.stamen.com/terrain/4/0/4.png
## Source : http://tile.stamen.com/terrain/4/1/4.png
## Source : http://tile.stamen.com/terrain/4/2/4.png
## Source : http://tile.stamen.com/terrain/4/0/5.png
## Source : http://tile.stamen.com/terrain/4/1/5.png
## Source : http://tile.stamen.com/terrain/4/2/5.png
ak_map_3857 <- ggmap_bbox_to_3857(ak_map)
- use inherit.aes = FALSE because there is stuff going on in the background
ggmap(ak_map_3857) +
geom_sf(data = pop_3857, aes(color = population), inherit.aes = F) +
scale_color_continuous(low = "khaki", high = "firebrick", labels = comma)
## Coordinate system already present. Adding new coordinate system, which will replace the existing one.
- there is a geom_raster that you can use if you have a raster basemap you want to use.
Visualize SF objects with leaflet
- pasted from book but not running or generating error message
epsg3338 <- leaflet::leafletCRS(
crsClass = "L.Proj.CRS",
code = "EPSG:3338",
proj4def = "+proj=aea +lat_1=55 +lat_2=65 +lat_0=50 +lon_0=-154 +x_0=0 +y_0=0 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs",
resolutions = 2^(16:7))
st_crs(pop_region_3338)
## Coordinate Reference System:
## User input: EPSG:3338
## wkt:
## PROJCRS["NAD83 / Alaska Albers",
## BASEGEOGCRS["NAD83",
## DATUM["North American Datum 1983",
## ELLIPSOID["GRS 1980",6378137,298.257222101,
## LENGTHUNIT["metre",1]]],
## PRIMEM["Greenwich",0,
## ANGLEUNIT["degree",0.0174532925199433]],
## ID["EPSG",4269]],
## CONVERSION["Alaska Albers (meters)",
## METHOD["Albers Equal Area",
## ID["EPSG",9822]],
## PARAMETER["Latitude of false origin",50,
## ANGLEUNIT["degree",0.0174532925199433],
## ID["EPSG",8821]],
## PARAMETER["Longitude of false origin",-154,
## ANGLEUNIT["degree",0.0174532925199433],
## ID["EPSG",8822]],
## PARAMETER["Latitude of 1st standard parallel",55,
## ANGLEUNIT["degree",0.0174532925199433],
## ID["EPSG",8823]],
## PARAMETER["Latitude of 2nd standard parallel",65,
## ANGLEUNIT["degree",0.0174532925199433],
## ID["EPSG",8824]],
## PARAMETER["Easting at false origin",0,
## LENGTHUNIT["metre",1],
## ID["EPSG",8826]],
## PARAMETER["Northing at false origin",0,
## LENGTHUNIT["metre",1],
## ID["EPSG",8827]]],
## CS[Cartesian,2],
## AXIS["easting (X)",east,
## ORDER[1],
## LENGTHUNIT["metre",1]],
## AXIS["northing (Y)",north,
## ORDER[2],
## LENGTHUNIT["metre",1]],
## USAGE[
## SCOPE["Topographic mapping (small scale)."],
## AREA["United States (USA) - Alaska."],
## BBOX[51.3,172.42,71.4,-129.99]],
## ID["EPSG",3338]]
pop_region_4326 <- pop_region_3338 %>% st_transform(crs = 4326)
m <- leaflet(options = leafletOptions(crs = epsg3338)) %>%
addPolygons(data = pop_region_4326,
fillColor = "gray",
weight = 1)
m
- We can add labels, legends, and a color scale.
pal <- colorNumeric(palette = "Reds", domain = pop_region_4326$total_pop)
m <- leaflet(options = leafletOptions(crs = epsg3338)) %>%
addPolygons(data = pop_region_4326,
fillColor = ~pal(total_pop),
weight = 1,
color = "black",
fillOpacity = 1,
label = ~region) %>%
addLegend(position = "bottomleft",
pal = pal,
values = range(pop_region_4326$total_pop),
title = "Total Population")
m
- We can also add the individual communities, with popup labels showing their population, on top of that!
pal <- colorNumeric(palette = "Reds", domain = pop_region_4326$total_pop)
m <- leaflet(options = leafletOptions(crs = epsg3338)) %>%
addPolygons(data = pop_region_4326,
fillColor = ~pal(total_pop),
weight = 1,
color = "black",
fillOpacity = 1) %>%
addCircleMarkers(data = pop_4326,
lat = ~lat,
lng = ~lng,
radius = ~log(population/500), # arbitrary scaling
fillColor = "gray",
fillOpacity = 1,
weight = 0.25,
color = "black",
label = ~paste0(pop_4326$city, ", population ", comma(pop_4326$population))) %>%
addLegend(position = "bottomleft",
pal = pal,
values = range(pop_region_4326$total_pop),
title = "Total Population")
m
- There is a lot more functionality to sf including the ability to intersect polygons, calculate distance, create a buffer, and more. Here are some more great resources and tutorials for a deeper dive into this great package: << See book for links! >>